
#include <limits.h>
#include <math.h>


namespace TC {

// Functions are defined in order of usage to assist inlining

inline void RandomMT::initialize( const uint32_t seed )
{
  // Initialize generator state with seed
  // See Knuth TAOCP Vol 2, 3rd Ed, p.106 for multiplier.
  // In previous versions, most significant bits (MSBs) of the seed affect
  // only MSBs of the state array.  Modified 9 Jan 2002 by Makoto Matsumoto.
  register uint32_t *s = m_state;
  register uint32_t *r = m_state;
  register int i = 1;
  *s++ = seed & 0xffffffffUL;
  for( ; i < N; ++i )
  {
    *s++ = ( 1812433253UL * ( *r ^ (*r >> 30) ) + i ) & 0xffffffffUL;
    r++;
  }
}

inline void RandomMT::reload()
{
  // Generate N new values in state
  // Made clearer and faster by Matthew Bellew (matthew.bellew@home.com)

  // in case enums are unsigned
  static TC_THREAD_LOCAL const int MmN = int(M) - int(N);  
  register uint32_t *p = m_state;
  register int i;
  for( i = N - M; i--; ++p )
    *p = twist( p[M], p[0], p[1] );
  for( i = M; --i; ++p )
    *p = twist( p[MmN], p[0], p[1] );
  *p = twist( p[MmN], p[0], m_state[0] );
  
  m_left = N, m_next = m_state;
}

inline void RandomMT::seed( const uint32_t oneSeed )
{
  // Seed the generator with a simple uint32_t
  initialize(oneSeed);
  reload();
}

inline void RandomMT::seed( uint32_t *const bigSeed, const uint32_t seedLength )
{
  // Seed the generator with an array of uint32_t's
  // There are 2^19937-1 possible initial states.  This function allows
  // all of those to be accessed by providing at least 19937 bits (with a
  // default seed length of N = 624 uint32_t's).  Any bits above the lower 32
  // in each element are discarded.
  // Just call seed() if you want to get array from /dev/urandom
  initialize(19650218UL);
  register int i = 1;
  register uint32_t j = 0;
  register int k = ( N > seedLength ? N : seedLength );
  for( ; k; --k )
  {
    m_state[i] =
    m_state[i] ^ ( (m_state[i-1] ^ (m_state[i-1] >> 30)) * 1664525UL );
    m_state[i] += ( bigSeed[j] & 0xffffffffUL ) + j;
    m_state[i] &= 0xffffffffUL;
    ++i;  ++j;
    if( i >= N ) { m_state[0] = m_state[N-1];  i = 1; }
    if( j >= seedLength ) j = 0;
  }
  for( k = N - 1; k; --k )
  {
    m_state[i] =
    m_state[i] ^ ( (m_state[i-1] ^ (m_state[i-1] >> 30)) * 1566083941UL );
    m_state[i] -= i;
    m_state[i] &= 0xffffffffUL;
    ++i;
    if( i >= N ) { m_state[0] = m_state[N-1];  i = 1; }
  }
  m_state[0] = 0x80000000UL;  // MSB is 1, assuring non-zero initial array
  reload();
}

inline void RandomMT::seed()
{
  // Seed the generator with an array from /dev/urandom if available
  // Otherwise use a hash of time() and clock() values
#if 0 
  // First try getting an array from /dev/urandom
  FILE* urandom = fopen( "/dev/urandom", "rb" );
  if( urandom )
  {
    uint32_t bigSeed[N];
    register uint32_t *s = bigSeed;
    register int i = N;
    register bool success = true;
    while( success && i-- )
      success = fread( s++, sizeof(uint32_t), 1, urandom );
    fclose(urandom);
    if( success ) { seed( bigSeed, N );  return; }
  }
#endif  
  // Was not successful, so use time() and clock() instead
  seed( Hash( time(NULL), clock() ) );
}

inline RandomMT::RandomMT( const uint32_t oneSeed )
  { seed(oneSeed); }

inline RandomMT::RandomMT( uint32_t *const bigSeed, const uint32_t seedLength )
  { seed(bigSeed,seedLength); }

inline RandomMT::RandomMT()
  { seed(); }

inline RandomMT::RandomMT( const RandomMT& o )
{
  register const uint32_t *t = o.m_state;
  register uint32_t *s = m_state;
  register int i = N;
  for( ; i--; *s++ = *t++ ) {}
  m_left = o.m_left;
  m_next = &m_state[N-m_left];
}

inline uint32_t RandomMT::randInt()
{
  // Pull a 32-bit integer from the generator state
  // Every other access function simply transforms the numbers extracted here
  
  if( m_left == 0 ) reload();
  --m_left;
  
  register uint32_t s1;
  s1 = *m_next++;
  s1 ^= (s1 >> 11);
  s1 ^= (s1 <<  7) & 0x9d2c5680UL;
  s1 ^= (s1 << 15) & 0xefc60000UL;
  return ( s1 ^ (s1 >> 18) );
}

inline uint32_t RandomMT::randInt( const uint32_t n )
{
  // Find which bits are used in n
  // Optimized by Magnus Jonsson (magnus@smartelectronix.com)
  uint32_t used = n;
  used |= used >> 1;
  used |= used >> 2;
  used |= used >> 4;
  used |= used >> 8;
  used |= used >> 16;
  
  // Draw numbers until one is found in [0,n]
  uint32_t i;
  do
    i = randInt() & used;  // toss unused bits to shorten search
  while( i > n );
  return i;
}

inline double RandomMT::rand()
  { return double(randInt()) * (1.0/4294967295.0); }

inline double RandomMT::rand( const double n )
  { return rand() * n; }

inline double RandomMT::randExc()
  { return double(randInt()) * (1.0/4294967296.0); }

inline double RandomMT::randExc( const double n )
  { return randExc() * n; }

inline double RandomMT::randDblExc()
  { return ( double(randInt()) + 0.5 ) * (1.0/4294967296.0); }

inline double RandomMT::randDblExc( const double n )
  { return randDblExc() * n; }

inline double RandomMT::rand53()
{
  uint32_t a = randInt() >> 5, b = randInt() >> 6;
  return ( a * 67108864.0 + b ) * (1.0/9007199254740992.0);  // by Isaku Wada
}

inline double RandomMT::randNorm( const double mean, const double stddev )
{
  // Return a real number from a normal (Gaussian) distribution with given
  // mean and standard deviation by polar form of Box-Muller transformation
  double x, y, r;
  do
  {
    x = 2.0 * rand() - 1.0;
    y = 2.0 * rand() - 1.0;
    r = x * x + y * y;
  }
  while ( r >= 1.0 || r == 0.0 );
  double s = sqrt( -2.0 * log(r) / r );
  return mean + x * s * stddev;
}

inline double RandomMT::operator()()
{
  return rand();
}

inline void RandomMT::save( uint32_t* saveArray ) const
{
  register const uint32_t *s = m_state;
  register uint32_t *sa = saveArray;
  register int i = N;
  for( ; i--; *sa++ = *s++ ) {}
  *sa = m_left;
}

inline void RandomMT::load( uint32_t *const loadArray )
{
  register uint32_t *s = m_state;
  register uint32_t *la = loadArray;
  register int i = N;
  for( ; i--; *s++ = *la++ ) {}
  m_left = *la;
  m_next = &m_state[N-m_left];
}


inline RandomMT& RandomMT::operator=( const RandomMT& o )
{
  if( this == &o ) return (*this);
  register const uint32_t *t = o.m_state;
  register uint32_t *s = m_state;
  register int i = N;
  for( ; i--; *s++ = *t++ ) {}
  m_left = o.m_left;
  m_next = &m_state[N-m_left];
  return (*this);
}


} // end namespace
